Using birth-death processes to infer tumor subpopulation structure from live-cell imaging drug screening data.

Tumor heterogeneity is a complex and widely recognized trait that poses significant challenges in developing effective cancer therapies. In particular, many tumors harbor a variety of subpopulations with distinct therapeutic response characteristics. Characterizing this heterogeneity by determining the subpopulation structure within a tumor enables more precise and successful treatment strategies. In our prior work, we developed PhenoPop, a computational framework for unravelling the drug-response subpopulation structure within a tumor from bulk high-throughput drug screening data. However, the deterministic nature of the underlying models driving PhenoPop restricts the model fit and the information it can extract from the data. As an advancement, we propose a stochastic model based on the linear birth-death process to address this limitation. Our model can formulate a dynamic variance along the horizon of the experiment so that the model uses more information from the data to provide a more robust estimation. In addition, the newly proposed model can be readily adapted to situations where the experimental data exhibits a positive time correlation. We test our model on simulated data (in silico) and experimental data (in vitro), which supports our argument about its advantages.

Statistical inference of the rates of cell proliferation and phenotypic switching in cancer.

Recent evidence suggests that nongenetic (epigenetic) mechanisms play an important role at all stages of cancer evolution. In many cancers, these mechanisms have been observed to induce dynamic switching between two or more cell states, which commonly show differential responses to drug treatments. To understand how these cancers evolve over time, and how they respond to treatment, we need to understand the state-dependent rates of cell proliferation and phenotypic switching. In this work, we propose a rigorous statistical framework for estimating these parameters, using data from commonly performed cell line experiments, where phenotypes are sorted and expanded in culture. The framework explicitly models the stochastic dynamics of cell division, cell death and phenotypic switching, and it provides likelihood-based confidence intervals for the model parameters. The input data can be either the fraction of cells or the number of cells in each state at one or more time points. Through a combination of theoretical analysis and numerical simulations, we show that when cell fraction data is used, the rates of switching may be the only parameters that can be estimated accurately. On the other hand, using cell number data enables accurate estimation of the net division rate for each phenotype, and it can even enable estimation of the state-dependent rates of cell division and cell death. We conclude by applying our framework to a publicly available dataset.

Limit theorems for the site frequency spectrum of neutral mutations in an exponentially growing population.

The site frequency spectrum (SFS) is a widely used summary statistic of genomic data, offering a simple means of inferring the evolutionary history of a population. Motivated by recent evidence for the role of neutral evolution in cancer, we examine the SFS of neutral mutations in an exponentially growing population. Whereas recent work has focused on the mean behavior of the SFS in this scenario, here, we investigate the first-order asymptotics of the underlying stochastic process. Using branching process techniques, we show that the SFS of a Galton-Watson process evaluated at a fixed time converges almost surely to a random limit. We also show that the SFS evaluated at the stochastic time at which the population first reaches a certain size converges in probability to a constant. Finally, we illustrate how our results can be used to construct consistent estimators for the extinction probability and the effective mutation rate of a birth-death process.

Statistical inference of the rates of cell proliferation and phenotypic switching in cancer.

Recent evidence suggests that nongenetic (epigenetic) mechanisms play an important role at all stages of cancer evolution. In many cancers, these mechanisms have been observed to induce dynamic switching between two or more cell states, which commonly show differential responses to drug treatments. To understand how these cancers evolve over time, and how they respond to treatment, we need to understand the state-dependent rates of cell proliferation and phenotypic switching. In this work, we propose a rigorous statistical framework for estimating these parameters, using data from commonly performed cell line experiments, where phenotypes are sorted and expanded in culture. The framework explicitly models the stochastic dynamics of cell division, cell death and phenotypic switching, and it provides likelihood-based confidence intervals for the model parameters. The input data can be either the fraction of cells or the number of cells in each state at one or more time points. Through a combination of theoretical analysis and numerical simulations, we show that when cell fraction data is used, the rates of switching may be the only parameters that can be estimated accurately. On the other hand, using cell number data enables accurate estimation of the net division rate for each phenotype, and it can even enable estimation of the state-dependent rates of cell division and cell death. We conclude by applying our framework to a publicly available dataset.

Phenotypic deconvolution in heterogeneous cancer cell populations using drug-screening data.

Tumor heterogeneity is an important driver of treatment failure in cancer since therapies often select for drug-tolerant or drug-resistant cellular subpopulations that drive tumor growth and recurrence. Profiling the drug-response heterogeneity of tumor samples using traditional genomic deconvolution methods has yielded limited results, due in part to the imperfect mapping between genomic variation and functional characteristics. Here, we leverage mechanistic population modeling to develop a statistical framework for profiling phenotypic heterogeneity from standard drug-screen data on bulk tumor samples. This method, called PhenoPop, reliably identifies tumor subpopulations exhibiting differential drug responses and estimates their drug sensitivities and frequencies within the bulk population. We apply PhenoPop to synthetically generated cell populations, mixed cell-line experiments, and multiple myeloma patient samples and demonstrate how it can provide individualized predictions of tumor growth under candidate therapies. This methodology can also be applied to deconvolution problems in a variety of biological settings beyond cancer drug response.

Exact site frequency spectra of neutrally evolving tumors: A transition between power laws reveals a signature of cell viability.

The site frequency spectrum (SFS) is a popular summary statistic of genomic data. While the SFS of a constant-sized population undergoing neutral mutations has been extensively studied in population genetics, the rapidly growing amount of cancer genomic data has attracted interest in the spectrum of an exponentially growing population. Recent theoretical results have generally dealt with special or limiting cases, such as considering only cells with an infinite line of descent, assuming deterministic tumor growth, or taking large-time or large-population limits. In this work, we derive exact expressions for the expected SFS of a cell population that evolves according to a stochastic branching process, first for cells with an infinite line of descent and then for the total population, evaluated either at a fixed time (fixed-time spectrum) or at the stochastic time at which the population reaches a certain size (fixed-size spectrum). We find that while the rate of mutation scales the SFS of the total population linearly, the rates of cell birth and cell death change the shape of the spectrum at the small-frequency end, inducing a transition between a 1/j2 power-law spectrum and a 1/j spectrum as cell viability decreases. We show that this insight can in principle be used to estimate the ratio between the rate of cell death and cell birth, as well as the mutation rate, using the site frequency spectrum alone. Although the discussion is framed in terms of tumor dynamics, our results apply to any exponentially growing population of individuals undergoing neutral mutations.

Understanding the role of phenotypic switching in cancer drug resistance.

The emergence of acquired drug resistance in cancer represents a major barrier to treatment success. While research has traditionally focused on genetic sources of resistance, recent findings suggest that cancer cells can acquire transient resistant phenotypes via epigenetic modifications and other non-genetic mechanisms. Although these resistant phenotypes are eventually relinquished by individual cells, they can temporarily 'save' the tumor from extinction and enable the emergence of more permanent resistance mechanisms. These observations have generated interest in the potential of epigenetic therapies for long-term tumor control or eradication. In this work, we develop a mathematical model to study how phenotypic switching at the single-cell level affects resistance evolution in cancer. We highlight unique features of non-genetic resistance, probe the evolutionary consequences of epigenetic drugs and explore potential therapeutic strategies. We find that even short-term epigenetic modifications and stochastic fluctuations in gene expression can drive long-term drug resistance in the absence of any bona fide resistance mechanisms. We also find that an epigenetic drug that slightly perturbs the average retention of the resistant phenotype can turn guaranteed treatment failure into guaranteed success. Lastly, we find that combining an epigenetic drug with an anti-cancer agent can significantly outperform monotherapy, and that treatment outcome is heavily affected by drug sequencing.

Glioblastoma Recurrence and the Role of O6-Methylguanine-DNA Methyltransferase Promoter Methylation.

Tumor recurrence in glioblastoma multiforme (GBM) is often attributed to acquired resistance to the standard chemotherapeutic agent, temozolomide (TMZ). Promoter methylation of the DNA repair gene MGMT (O6-methylguanine-DNA methyltransferase) has been associated with sensitivity to TMZ, whereas increased expression of MGMT has been associated with TMZ resistance. Clinical studies have observed a downward shift in MGMT methylation percentage from primary to recurrent stage tumors; however, the evolutionary processes that drive this shift and more generally the emergence and growth of TMZ-resistant tumor subpopulations are still poorly understood. Here, we develop a mathematical model, parameterized using clinical and experimental data, to investigate the role of MGMT methylation in TMZ resistance during the standard treatment regimen for GBM-surgery, chemotherapy, and radiation. We first found that the observed downward shift in MGMT promoter methylation status between detection and recurrence cannot be explained solely by evolutionary selection. Next, our model suggests that TMZ has an inhibitory effect on maintenance methylation of MGMT after cell division. Finally, incorporating this inhibitory effect, we study the optimal number of TMZ doses per adjuvant cycle for patients with GBM with high and low levels of MGMT methylation at diagnosis.

A Robust Optimization Approach to Cancer Treatment under Toxicity Uncertainty.

The design of optimal protocols plays an important role in cancer treatment. However, in clinical applications, the outcomes under the optimal protocols are sensitive to variations of parameter settings such as drug effects and the attributes of age, weight, and health conditions in human subjects. One approach to overcoming this challenge is to formulate the problem of finding an optimal treatment protocol as a robust optimization problem (ROP) that takes parameter uncertainty into account. In this chapter, we describe a method to model toxicity uncertainty. We then apply a mixed integer ROP to derive the optimal protocols that minimize the cumulative tumor size. While our method may be applied to other cancers, in this work we focus on the treatment of chronic myeloid leukemia (CML) with tyrosine kinase inhibitors (TKI). For simplicity, we focus on one particular mode of toxicity arising from TKI therapy, low blood cell counts, in particular low absolute neutrophil count (ANC). We develop optimization methods for locating optimal treatment protocols assuming that the rate of decrease of ANC varies within a given interval. We further investigated the relationship between parameter uncertainty and optimal protocols. Our results suggest that the dosing schedule can significantly reduce tumor size without recurrence in 360 weeks while insuring that toxicity constraints are satisfied for all realizations of uncertain parameters.

Combination therapeutics of Nilotinib and radiation in acute lymphoblastic leukemia as an effective method against drug-resistance.

Philadelphia chromosome-positive (Ph+) acute lymphoblastic leukemia (ALL) is characterized by a very poor prognosis and a high likelihood of acquired chemo-resistance. Although tyrosine kinase inhibitor (TKI) therapy has improved clinical outcome, most ALL patients relapse following treatment with TKI due to the development of resistance. We developed an in vitro model of Nilotinib-resistant Ph+ leukemia cells to investigate whether low dose radiation (LDR) in combination with TKI therapy overcome chemo-resistance. Additionally, we developed a mathematical model, parameterized by cell viability experiments under Nilotinib treatment and LDR, to explain the cellular response to combination therapy. The addition of LDR significantly reduced drug resistance both in vitro and in computational model. Decreased expression level of phosphorylated AKT suggests that the combination treatment plays an important role in overcoming resistance through the AKT pathway. Model-predicted cellular responses to the combined therapy provide good agreement with experimental results. Augmentation of LDR and Nilotinib therapy seems to be beneficial to control Ph+ leukemia resistance and the quantitative model can determine optimal dosing schedule to enhance the effectiveness of the combination therapy.

Quantifying the Dynamics of Field Cancerization in Tobacco-Related Head and Neck Cancer: A Multiscale Modeling Approach.

High rates of local recurrence in tobacco-related head and neck squamous cell carcinoma (HNSCC) are commonly attributed to unresected fields of precancerous tissue. Because they are not easily detectable at the time of surgery without additional biopsies, there is a need for noninvasive methods to predict the extent and dynamics of these fields. Here, we developed a spatial stochastic model of tobacco-related HNSCC at the tissue level and calibrated the model using a Bayesian framework and population-level incidence data from the Surveillance, Epidemiology, and End Results (SEER) registry. Probabilistic model analyses were performed to predict the field geometry at time of diagnosis, and model predictions of age-specific recurrence risks were tested against outcome data from SEER. The calibrated models predicted a strong dependence of the local field size on age at diagnosis, with a doubling of the expected field diameter between ages at diagnosis of 50 and 90 years, respectively. Similarly, the probability of harboring multiple, clonally unrelated fields at the time of diagnosis was found to increase substantially with patient age. On the basis of these findings, we hypothesized a higher recurrence risk in older than in younger patients when treated by surgery alone; we successfully tested this hypothesis using age-stratified outcome data. Further clinical studies are needed to validate the model predictions in a patient-specific setting. This work highlights the importance of spatial structure in models of epithelial carcinogenesis and suggests that patient age at diagnosis may be a critical predictor of the size and multiplicity of precancerous lesions. Cancer Res; 76(24); 7078-88. ©2016 AACR.

Optimized Treatment Schedules for Chronic Myeloid Leukemia.

Over the past decade, several targeted therapies (e.g. imatinib, dasatinib, nilotinib) have been developed to treat Chronic Myeloid Leukemia (CML). Despite an initial response to therapy, drug resistance remains a problem for some CML patients. Recent studies have shown that resistance mutations that preexist treatment can be detected in a substantial number of patients, and that this may be associated with eventual treatment failure. One proposed method to extend treatment efficacy is to use a combination of multiple targeted therapies. However, the design of such combination therapies (timing, sequence, etc.) remains an open challenge. In this work we mathematically model the dynamics of CML response to combination therapy and analyze the impact of combination treatment schedules on treatment efficacy in patients with preexisting resistance. We then propose an optimization problem to find the best schedule of multiple therapies based on the evolution of CML according to our ordinary differential equation model. This resulting optimization problem is nontrivial due to the presence of ordinary different equation constraints and integer variables. Our model also incorporates drug toxicity constraints by tracking the dynamics of patient neutrophil counts in response to therapy. We determine optimal combination strategies that maximize time until treatment failure on hypothetical patients, using parameters estimated from clinical data in the literature.

Optimal treatment and stochastic modeling of heterogeneous tumors.

UNLABELLED: In this work we review past articles that have mathematically studied cancer heterogeneity and the impact of this heterogeneity on the structure of optimal therapy. We look at past works on modeling how heterogeneous tumors respond to radiotherapy, and take a particularly close look at how the optimal radiotherapy schedule is modified by the presence of heterogeneity. In addition, we review past works on the study of optimal chemotherapy when dealing with heterogeneous tumors. REVIEWERS: This article was reviewed by Thomas McDonald, David Axelrod, and Leonid Hanin.

A mathematical model of tumor growth and its response to single irradiation.

BACKGROUND: Mathematical modeling of biological processes is widely used to enhance quantitative understanding of bio-medical phenomena. This quantitative knowledge can be applied in both clinical and experimental settings. Recently, many investigators began studying mathematical models of tumor response to radiation therapy. We developed a simple mathematical model to simulate the growth of tumor volume and its response to a single fraction of high dose irradiation. The modelling study may provide clinicians important insights on radiation therapy strategies through identification of biological factors significantly influencing the treatment effectiveness. METHODS: We made several key assumptions of the model. Tumor volume is composed of proliferating (or dividing) cancer cells and non-dividing (or dead) cells. Tumor growth rate (or tumor volume doubling time) is proportional to the ratio of the volumes of tumor vasculature and the tumor. The vascular volume grows slower than the tumor by introducing the vascular growth retardation factor, θ. Upon irradiation, the proliferating cells gradually die over a fixed time period after irradiation. Dead cells are cleared away with cell clearance time. The model was applied to simulate pre-treatment growth and post-treatment radiation response of rat rhabdomyosarcoma tumors and metastatic brain tumors of five patients who were treated with Gamma Knife stereotactic radiosurgery (GKSRS). RESULTS: By selecting appropriate model parameters, we showed the temporal variation of the tumors for both the rat experiment and the clinical GKSRS cases could be easily replicated by the simple model. Additionally, the application of our model to the GKSRS cases showed that the α-value, which is an indicator of radiation sensitivity in the LQ model, and the value of θ could be predictors of the post-treatment volume change. CONCLUSIONS: The proposed model was successful in representing both the animal experimental data and the clinically observed tumor volume changes. We showed that the model can be used to find the potential biological parameters, which may be able to predict the treatment outcome. However, there is a large statistical uncertainty of the result due to the small sample size. Therefore, a future clinical study with a larger number of patients is needed to confirm the finding.

Optimal radiotherapy dose schedules under parametric uncertainty.

We consider the effects of parameter uncertainty on the optimal radiation schedule in the context of the linear-quadratic model. Our interest arises from the observation that if inter-patient variability in normal and tumor tissue radiosensitivity or sparing factor of the organs-at-risk (OAR) are not accounted for during radiation scheduling, the performance of the therapy may be strongly degraded or the OAR may receive a substantially larger dose than the allowable threshold. This paper proposes a stochastic radiation scheduling concept to incorporate inter-patient variability into the scheduling optimization problem. Our method is based on a probabilistic approach, where the model parameters are given by a set of random variables. Our probabilistic formulation ensures that our constraints are satisfied with a given probability, and that our objective function achieves a desired level with a stated probability. We used a variable transformation to reduce the resulting optimization problem to two dimensions. We showed that the optimal solution lies on the boundary of the feasible region and we implemented a branch and bound algorithm to find the global optimal solution. We demonstrated how the configuration of optimal schedules in the presence of uncertainty compares to optimal schedules in the absence of uncertainty (conventional schedule). We observed that in order to protect against the possibility of the model parameters falling into a region where the conventional schedule is no longer feasible, it is required to avoid extremal solutions, i.e. a single large dose or very large total dose delivered over a long period. Finally, we performed numerical experiments in the setting of head and neck tumors including several normal tissues to reveal the effect of parameter uncertainty on optimal schedules and to evaluate the sensitivity of the solutions to the choice of key model parameters.

Minimizing metastatic risk in radiotherapy fractionation schedules.

Metastasis is the process by which cells from a primary tumor disperse and form new tumors at distant anatomical locations. The treatment and prevention of metastatic cancer remains an extremely challenging problem. This work introduces a novel biologically motivated objective function to the radiation optimization community that takes into account metastatic risk instead of the status of the primary tumor. In this work, we consider the problem of developing fractionated irradiation schedules that minimize production of metastatic cancer cells while keeping normal tissue damage below an acceptable level. A dynamic programming framework is utilized to determine the optimal fractionation scheme. We evaluated our approach on a breast cancer case using the heart and the lung as organs-at-risk (OAR). For small tumor [Formula: see text] values, hypo-fractionated schedules were optimal, which is consistent with standard models. However, for relatively larger [Formula: see text] values, we found the type of schedule depended on various parameters such as the time when metastatic risk was evaluated, the [Formula: see text] values of the OARs, and the normal tissue sparing factors. Interestingly, in contrast to standard models, hypo-fractionated and semi-hypo-fractionated schedules (large initial doses with doses tapering off with time) were suggested even with large tumor α/ß values. Numerical results indicate the potential for significant reduction in metastatic risk.

An Evolutionary Approach for Identifying Driver Mutations in Colorectal Cancer.

The traditional view of cancer as a genetic disease that can successfully be treated with drugs targeting mutant onco-proteins has motivated whole-genome sequencing efforts in many human cancer types. However, only a subset of mutations found within the genomic landscape of cancer is likely to provide a fitness advantage to the cell. Distinguishing such "driver" mutations from innocuous "passenger" events is critical for prioritizing the validation of candidate mutations in disease-relevant models. We design a novel statistical index, called the Hitchhiking Index, which reflects the probability that any observed candidate gene is a passenger alteration, given the frequency of alterations in a cross-sectional cancer sample set, and apply it to a mutational data set in colorectal cancer. Our methodology is based upon a population dynamics model of mutation accumulation and selection in colorectal tissue prior to cancer initiation as well as during tumorigenesis. This methodology can be used to aid in the prioritization of candidate mutations for functional validation and contributes to the process of drug discovery.

The Impact of Microenvironmental Heterogeneity on the Evolution of Drug Resistance in Cancer Cells.

Therapeutic resistance arises as a result of evolutionary processes driven by dynamic feedback between a heterogeneous cell population and environmental selective pressures. Previous studies have suggested that mutations conferring resistance to epidermal growth factor receptor (EGFR) tyrosine kinase inhibitors (TKI) in non-small-cell lung cancer (NSCLC) cells lower the fitness of resistant cells relative to drug-sensitive cells in a drug-free environment. Here, we hypothesize that the local tumor microenvironment could influence the magnitude and directionality of the selective effect, both in the presence and absence of a drug. Using a combined experimental and computational approach, we developed a mathematical model of preexisting drug resistance describing multiple cellular compartments, each representing a specific tumor environmental niche. This model was parameterized using a novel experimental dataset derived from the HCC827 erlotinib-sensitive and -resistant NSCLC cell lines. We found that, in contrast to in the drug-free environment, resistant cells may hold a fitness advantage compared to parental cells in microenvironments deficient in oxygen and nutrients. We then utilized the model to predict the impact of drug and nutrient gradients on tumor composition and recurrence times, demonstrating that these endpoints are strongly dependent on the microenvironment. Our interdisciplinary approach provides a model system to quantitatively investigate the impact of microenvironmental effects on the evolutionary dynamics of tumor cells.

Multifocality and recurrence risk: a quantitative model of field cancerization.

Primary tumors often emerge within genetically altered fields of premalignant cells that appear histologically normal but have a high chance of progression to malignancy. Clinical observations have suggested that these premalignant fields pose high risks for emergence of recurrent tumors if left behind after surgical removal of the primary tumor. In this work, we develop a spatio-temporal stochastic model of epithelial carcinogenesis, combining cellular dynamics with a general framework for multi-stage genetic progression to cancer. Using the model, we investigate how various properties of the premalignant fields depend on microscopic cellular properties of the tissue. In particular, we provide analytic results for the size-distribution of the histologically undetectable premalignant fields at the time of diagnosis, and investigate how the extent and the geometry of these fields depend upon key groups of parameters associated with the tissue and genetic pathways. We also derive analytical results for the relative risks of local vs. distant secondary tumors for different parameter regimes, a critical aspect for the optimal choice of post-operative therapy in carcinoma patients. This study contributes to a growing literature seeking to obtain a quantitative understanding of the spatial dynamics in cancer initiation.

Mathematical modeling of PDGF-driven glioblastoma reveals optimized radiation dosing schedules.

Glioblastomas (GBMs) are the most common and malignant primary brain tumors and are aggressively treated with surgery, chemotherapy, and radiotherapy. Despite this treatment, recurrence is inevitable and survival has improved minimally over the last 50 years. Recent studies have suggested that GBMs exhibit both heterogeneity and instability of differentiation states and varying sensitivities of these states to radiation. Here, we employed an iterative combined theoretical and experimental strategy that takes into account tumor cellular heterogeneity and dynamically acquired radioresistance to predict the effectiveness of different radiation schedules. Using this model, we identified two delivery schedules predicted to significantly improve efficacy by taking advantage of the dynamic instability of radioresistance. These schedules led to superior survival in mice. Our interdisciplinary approach may also be applicable to other human cancer types treated with radiotherapy and, hence, may lay the foundation for significantly increasing the effectiveness of a mainstay of oncologic therapy. PAPERCLIP: